home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Programmer Power Tools
/
Programmer Power Tools.iso
/
progjrn
/
pj_6_6.arc
/
FCODE.ARC
/
CHILC.FOR
next >
Wrap
Text File
|
1987-08-31
|
41KB
|
576 lines
C./ ADD NAME=CHILC
C PROGRAM TO CALCULATE, PRINT, AND PLOT TNE MAGNETIC
C SUSCEPTIBILITY OF HSP METALS.
C
IMPLICIT REAL*8 (A-H,P-Z)
REAL FUN(200,5),NAM(1000),NBOUND(10,2)
DIMENSION S(1300,3),TM(2,2600),T(2),C(4),S1(3),TE(2,4),TO(2,5)
COMMON/B/ CRIT, DE ,NTI
COMMON/A/NA, NS, KL, ISIG
COMMON/C/EPS, DLE
COMMON/D/ S12, S11
REAL*8 KL,KAPPA,KAPRM
C
PI=3.1415927D0
open (21,FILE='FOR21.DAT')
open (5,FILE='CHILC.OUT')
7 READ (21,4100) KODE
400 FORMAT(' ENTER RESTART KODE ',I5,/)
print 400, KODE
IF (KODE.LT.1) KODE = 1
IF (KODE.GT.6) STOP
GO TO (2, 1, 5, 6, 22, 35), KODE
2 READ (21,600) NTI, CRIT, DT, CN, I, NI, NR
500 FORMAT (' ENTER NTI, CRIT, DT, CN, I, NI, NR',I5,3F10.5,3I5,/)
print 500, NTI, CRIT, DT, CN, I, NI, NR
600 FORMAT (I5, 3F10.5, 3I5 )
C NTI IS NUMBER OF SIMPSON INTERVAL IN INTEGRATION FOR TM.
C CRIT IS WIDTH (IN KT) OF FERMI FUNCTION USED.
C DT IS RANGE (IN KT) ABOVE FERMI LEVEL WHERE TM CALCULATED.
C CN CONTROLS SUBDIVIDING OF SIMPSON INTERVALS IN E INTEGRATION.
C IF CN = .0, ENTERED VALUES OF I, NI, NR ARE USEL.
C I AND NI CONTROL RRNGE OF SUBDIVISION, FROM L=I TO L=I+2*NI.
C NR IS THE NUMWER OF SUBDIVISIONS PER INTERVAL, SDE=DE/(2NR).
C IF NTI = 0, STANDARD VALUES OF ABOVE ARE ENTERED, SEE BELOW.
IF (NTI.LT.0) GO TO 7
IF (NTI.GT.0) GO TO 1
NTI = 20
CRIT = 10.
DT = 10.
CN = .0
I = 0
NI = 0
NR = 1
C
C SET UP S VECTORS
1 READ (21,2000) NE, DELTA, EC, BET, ISIG
1000 FORMAT (' ENTER NE, DELTA, EC, BETA, ISIG ',I5,3F10.5,I5,/)
print 1000, NE, DELTA, EC, BET, ISIG
C DELTA, EC, AND BET=BETA=B/A ARE BENNETT-FALICOV BAND PARAMETERS.
C NE IS THE NUMBER OF SIMPSON INTERVALS, DE = (EMAX-EMIN)/(2NE).
C E2/(2DE) AND E3/(2DE) MUST BE INTEGERS, LATER IN THE
C PROGRAM, NE BECOMES 2NE+1.
C IN THIS SECTION IF ISIG GT 1, GET TRACE AND S1,S2,S;
2000 FORMAT (I5, 3F10.5, I5)
IF (NE.LE.0) GO TO 2
B2 = BET**.66666666D0
B4 = B2**2
RP3 = (B4 + 2./B2)/3.D0
BE2 = BET**2
E3 = EC + DELTA/3.D0
E2 = EC - DELTA/3.D0
EMIN = .0
IF (E2.LT..0) EMIN = E2
EMAX = .0
IF (E3.GT.0.) EMAX = E3
EMID = E2
IF (EMID.EQ.EMIN) EMID = .0
IF (EMID.EQ.EMAX) EMID = E3
EL = .0
IF (EC.EQ.0.) GO TO 4
EL=8.*EC**3/((((4.*RP3/(3.*B4))-1.)*EC**2+DELTA**2/9.D0)*27.*BE2)
IF (BET.EQ.1.) GO TO 4
21 KAPPA = EL*(EL-EC-DELTA/3.)*(EL-EC+DELTA/3.)
KAPRM = 2.*EL*(EL-EC) + ((EL-EC)**2-(DELTA/3.)**2)
H3 = RP3*EL - 2.*EC/(3.*B2)
DEL = -(KAPPA-H3**3)/(KAPRM-3.*RP3*H3**2)
EL = EL + DEL
IF (DABS(DEL/(EL-EMID)).GT..001) GO TO 21
4 print 1500, EL, E2, E3
1500 FORMAT (' EL =',E12.4,' E2 = ',E12.4,' E3 = ',E12.4,/ )
DE = (EMAX - EMIN)/(2.*NE)
NE = 2*NE + 1
9 IF (NE.GT.1300) GO TO 90
IT=0
NIT=0
NB=0
IF (CN.EQ.0.) GO TO 101
C PICK SUBDIVISION PARAMETERS
I=0
NI=0
NR=1
101 C(1)=(EMAX-EMID)/DE
C(2)=(EMID-EMIN)/DE
C(3)=(EL-EMID)/DE
TEMP=DABS(C(3))
IF (C(1).GT.C(2)) GOTO 104
IF (C(1).LT.1.) GOTO 11
IF (C(1).GE.CN) GOTO 104
NR = CN/C(1)+1.
NI = C(1)+.1
I = NE-2*NI
IF (I.GE.1) GOTO 107
NI = NI+(I-1)/2
I = 1
GOTO 107
11 IF (TEMP.LE.2.) GOTO 95
EMAX=EMAX - 2.*DE
NE = NE-2
IF (CN.EQ.0.) GOTO 109
C(1) = C(1)-2.
C(3) = C(3)+2.
TEMP = TEMP-2.
I = NE-2
NI=TEMP/4.
IF (NI.LT.2) NI=2
NR = CN
GOTO 107
104 IF (C(2).LT.1.) GOTO 12
IF (C(2).GT.CN) GOTO 107
NR = CN/C(2)+1.
NI = C(2)+.1
I = 1
GOTO 107
12 IF (TEMP.LE.2.) GOTO 95
EMIN = EMIN+2.*DE
NE = NE-2
IF (CN.EQ.0.) GOTO 109
C(2) = C(2)-2.
C(3) = C(3)-2.
TEMP = TEMP-2.
I = 1
NI = TEMP/4.
IF (NI.LT.2) NI=2
NR = CN
107 IF (TEMP.GE.CN.OR.CN.EQ.0.) GOTO 109
IF (TEMP.LT..5/CN) GOTO 109
NRT = CN/TEMP+.999
NIT = IFIX(SNGL(TEMP+1.))+1
IT = -1-2*(IFIX(SNGL(TEMP))/4)
IF (C(3).LT.0.) IT = 2-2*NIT-IT
IT = IFIX(SNGL(C(2)+.1))+IT
IF (NRT.GT.NR) NR = NRT
IF (IT.LT.1) IT = 1
IF (I.EQ.0) I = IT
IF (I.LE.IT) GOTO 102
NI = NI+(I-IT)/2
I = IT
102 IF (I+2*NI.LT.IT+2*NIT) NI = (IT-I)/2+NIT
109 print 1501, I, NI, NR, NB
C
C ENERGY INTEGRATION
E = EMIN
DO 111 L = 1,NE
DO 111 K = 1,3
111 S(L,K)=.0
IF (ISIG.GT.1) print 1501, I, NI, NR, NB
1501 FORMAT(' I, NI, NR, NB = ',4I4,/)
NB = 1
IF (I.EQ.1) NB = NR
IF (ISIG.GT.1) print 1501, I, NI, NR, NB
SDE = DE/NB
DLE = .0
CALL SS( E, SDE, S1, DELTA, EC, BET)
NEM1= NE-1
DO 3 J = 2, NEM1, 2
NB = 1
IF (J.GT.I.AND.J.LT.I+2*NI) NB = NR
IF (ISIG.GT.1) print 1502, I, NI, NR, NB, J
1502 FORMAT (5I4)
SDE = DE/NB
DO 112 K = 1,3
112 S(J-1,K) = S(J-1,K)+S1(K)/(2.*NB)
IF (ISIG.GT.1) print 6000, E, ( S(J-1,K),K = 1,3)
NOSC = 1
A = -1.+1./NB
NB2 = NB*2
DO 115 M = 1,NB2
E = E+SDE
C IF (J.EQ.NE-1.AND.M.EQ.2*NB) E = EMAX
C IF (DABS(E-EMID).LT..5*SDE.AND.M.EQ.2*NB) E = EMID
C DLE = .0
EPS = E-EL
IF (EPS.GT.-.5*SDE.AND.EPS.LE..5*SDE) DLE = SDE
IF (DLE.EQ.0.D0) GOTO 114
DO 113 K = 1,3
113 TM(1,K) = S1(K)
IF (DLE.NE.0.D0.AND.DABS(E-EMID).LT..5*SDE) GOTO 17
114 CALL SS( E, SDE, S1, DELTA, EC, BET)
IF (DLE.EQ.0.D0) GOTO 15
C
C CALCULATE CONTRIBUTIONNOF SINGULARITY AT E = EL
DO 13 K = 1,3
13 C(K) = S1(K)
DLE = .0D0
CALL SS( E+SDE, SDE, S1, DELTA, EC, BET)
IF (DABS(E-EMID).GT..1*SDE) GOTO 119
S1(1) = S1(1)+S11
S1(2) = S1(2)+S12
119 TEMP = EPS/SDE
DO 14 K = 1,3
IF (MOD(M,2).NE.0) S1(K)=(2.*S1(K)+2.*TM(1,K)+(3.D0*TEMP
1 *DLOG((1.+TEMP)/(1.-TEMP))-6.D0)*C(K))/4.D0
IF (MOD(M,2).EQ.0) S1(K) = (S1(K)+TM(1,K)+(3.D0*TEMP
1 *DLOG((1.+TEMP)/(1.-TEMP))+.5D0*DLOG((4.-TEMP**2)/(1.-TEMP**2)
2 )-6.D0)*C(K))/2.
14 CONTINUE
GOTO 15
C
C CALCULATE CONTRIBUTION OF SINGULARITY AT E = EL = EMID
17 DLE = .0D0
DTE = DABS(EMID-2.*EC/(3.*RP3*B2))
IF (DTE.GT..5*SDE) DTE = .5*SDE
IF (DTE.LT..05*SDE) DTE = .05*SDE
TW3 = 2.**.33333333D0
TEMP = 3.D0/(2.D0*TW3)-1.
TEMP3 = (SDE/DTE)**.33333333D0-2.+DTE/SDE+(1.-DTE/SDE)/TW3
TEMP1 = 3.D0*TEMP/TEMP3
TEMP2 = (1.+3.D0*TEMP*(DTE/SDE-2.)/TEMP3)
TEMP3 = -.5D0+3.*TEMP*(1.-DTE/SDE)/TEMP3
CALL SS ( E-2.*SDE, SDE, C, DELTA, EC, BET)
DO 117 K = 1,3
117 TM(1,K) = TEMP2*S1(K)+TEMP3*C(K)
CALL SS ( E+SDE, SDE, S1, DELTA, EC, BET)
CALL SS ( E+2.*SDE, SDE, C, DELTA, EC, BET)
DO 121 K = 1,3
121 TM(1,K) = TM(1,K)+TEMP2*S1(K)+TEMP3*C(K)
CALL SS ( E-DTE, DTE, S1, DELTA, EC, BET)
CALL SS ( E+DTE, DTE, C, DELTA, EC, BET)
DO 118 K = 1,3
118 S1(K) = TM(1,K)+TEMP1*(S1(K)+C(K))
C
15 TEMP = (3.D0+NOSC)/(4.D0*NB)
IF ( M.EQ.2*NB) TEMP = TEMP/2.
IF (ISIG.GT.1) print 1503, M, NOSC, A, TEMP
1503 FORMAT (' M, NOSC, A, TEMP ', 2I4, 2E12.4)
DO 116 K = 1,3
S(J-1,K) = S(J-1,K)-TEMP*A*(1.-A)*S1(K)
S(J,K) = S(J,K)+TEMP*(1.-A**2)*S1(K)
116 S(J+1,K) = S(J+1,K)+TEMP*A*(1.+A)*S1(K)
IF (ISIG.GT.1) print 6000, E, ( S(J-1,K), K = 1,3)
IF (ISIG.GT.1) print 6000, E, ( S(J,K), K = 1,3)
IF (ISIG.GT.1) print 6000, E, ( S(J+1,K), K = 1,3)
A = A+1./NB
115 NOSC = -NOSC
IF (DABS(E-EMID).GT..5*SDE.OR.DABS(EPS).LT..5*SDE) GOTO 3
C
C STEP FUNCTION CORRECTOR
TEMP = S11*DSIGN(1.D0,EL-EMID)
S(J+1,1) = S(J+1,1)-TEMP/(2.*NB)
S1(1) = S1(1)+TEMP
TEMP = S12*DSIGN(1.D0,EL-EMID)
S(J+1,2) = S(J+1,2)-TEMP/(2.*NB)
S1(2) = S1(2)+TEMP
C
3 CONTINUE
DO 120 K = 1,3
S(1,K) = 2.*S(1,K)
120 S(NE,K) = 2.*S(NE,K)
C
C PLOT OR PRINT S; ISIG = 0, PLOT S(M); ISIG > 0. PRINT;
C ISIG < 0. PRINT AND PLOT.
READ (21,4100) ISIG,M
3000 FORMAT (' ENTER ISIG,M ',2I5,/)
print 3000, ISIG,M
4100 FORMAT ( 2I5 )
IF (ISIG) 5,5,6
5 IF (M.LT.1.OR.M.GT.3) GOTO 20
SMIN = S(1,M)
SMAX = SMIN
DO 16 L = 2,NE
SMIN = DMIN1(SMIN,S(L,M))
16 SMAX = DMAX1(SMAX,S(L,M))
print 4200, SMIN, SMAX
4200 FORMAT (' SMIN =',E12.4,' SMAX =',E12.4/' ENTER SMIN,SMAX '/)
C READ (21,4000), SMIN, SMAX
4000 FORMAT (2F10.5)
C CALL PLOT (EMIN,.0,12)
C CALL PLOT (.0,SMIN,13)
C CALL PLOT (.0,SMIN,12)
C CALL PLOT (.0,SMAX,12)
C CALL TTY
6 IF (ISIG.EQ.3) print 5000
5000 FORMAT(' ',5X,'E',8X,'S(1)',8X,'S(2)',8X,'N(E)')
E = EMIN
DO 10 L = 1,NE
IF (ISIG.EQ.3) print 6000, E, ( S(L,K), K = 1,3)
6000 FORMAT (F12.6, 3E12.4)
C IF (ISIG.LE.0.AND.L.EQ.1) CALL PLOT(EMIN,S(L,M),13)
C IF (ISIG.LE.0) CALL PLOT ( E, S(L,M), 12)
C IF (ISIG.LT.0) CALL TTY
10 E = E+DE
C CALL TTY
IF (ISIG.LE.0) GOTO 5
C
C PREPARE S WITH SIMPSON FACTORS
20 OSC = DE/3.D0
DO 8 L = 1,NE
OSC = -OSC
SIMP = DE+OSC
IF (L.EQ.1.OR.L.EQ.NE) SIMP = DE/3.D0
DO 8 K = 1,3
8 S(L,K) = S(L,K)*SIMP
C
LGR = 1
C SET UP TM'S
22 READ (21,8000) EFMIN, EFMAX, TRY, A, B
7000 FORMAT (' ENTER EFMIN, EFMAX, TRY, ABAR, B ',5F10.6,/)
print 7000, EFMIN, EFMAX, TRY, A, B
8000 FORMAT (5F10.6)
C EFMIN, EFMAX ARE MINIMUM AND MAXIMUM FERMI LEVELS
C TRY IS THE TEMPERATURE (IN ENERGY UNITS)
C ABAR AND B ARE BAND PARAMETERS (NOT THE B/A =BETA)
IF (B.EQ.0.) GOTO 1
C FMUL = 1.71732D-6*A**2/DSQRT(B)
FMUL = 1.71732D-6*A**2/DSQRT(B) /7.14
FNE = 1.8769D4*RP3/A**4
NM = -(EFMIN-EMIN)/DE
EFMIN = EMIN-NM*DE
NM = (EFMAX-EMIN)/DE
EFMAX = EMIN+NM*DE
NF = (EFMAX-EFMIN)/DE+1.1
NT = (EFMAX-EMIN+DT*TRY)/DE+2.1
IF (NT.GT.NE+NF) NT = NE+NF
IF (NT.GT.2600) GOTO 93
IF (NT.LT.1) GOTO 35
EMU = EMIN-EFMAX
DO 30 L = 1,NT
CALL TN( 1.D0, TRY, EMU, T)
IF (ISIG.EQ.4) print 6000, EMU, ( T(K), K = 1,2)
DO 25 K = 1,2
25 TM(K,L) = T(K)*FMUL
30 EMU = EMU+DE
IF (TRY.NE.0.) GOTO 35
CALL TONE(2.D0,TO)
SM = DSQRT(DE)*FMUL
DO 32 K = 1,2
DO 31 N = 1,4
31 TO(K,N) = TO(K,N)*SM
IF (ISIG.GT.2) print 6000, SM, (TO(K,L), L=1,4)
32 SM = SM/DE
OSC = 1.D0/3.D0
DO 34 N = 1,4
DO 33 K = 1,2
TE(K,N) = TO(K,N)/(1.-OSC)
33 TO(K,N) = TO(K,N)/(1.+OSC)
34 OSC = -OSC
DO 28 K = 1,2
TO(K,5) = TO(K,4)
TO(K,4) = TO(K,3)
TO(K,3) = TO(K,2)
TO(K,2) = TO(K,1)+3.D0*TM(K,NT-3)/8.D0
TO(K,1) = 5.D0*TM(K,NT-4)/4.D0
28 IF (ISIG.GT.2) print 6000, OSC, (TO(K,L), L = 1,5)
DO 29 K = 1,2
TM(K,NT-3) = TE(K,1)+TM(K,NT-3)/2.
TM(K,NT-2) = TE(K,2)
TM(K,NT-1) = TE(K,3)
29 TM(K,NT) = TE(K,4)
C
C PLOT OR PRINT C'S AND CHI = C(3), N(E) = C(4)
C IF ISIG = 0, PLOT C(M)
C IF ISIG > 0, PRINT ALL C'S
C IF ISIG < 0, PRINT ALL C'S AND PLOT C(M)
35 READ (21,8800) ISIG, M,NPR
8500 FORMAT(' ENTER ISIG, M, NPR',3I5,/)
print 8500, ISIG, M, NPR
8800 FORMAT (3I5)
IF (M.LE.0.OR.M.GT.4) GOTO 22
IF (ISIG) 36,36,37
36 print 9000
9000 FORMAT (' ENTER CMIN, CMAX '/)
READ (21,4000) CMIN, CMAX
DO 52 K=1,5
NBOUND(K,1) = CMIN
52 NBOUND(K,2) = CMAX
9700 FORMAT(9X,'EF',4X,'C(1)',7X,'C(2)',7X,'CHI',8X,'N(E)',/)
37 IF (ISIG.GE.1) print 9700
C
C COMRUTE C'S
EF = EFMIN
DO 50 N = 1,NF
DO 38 K = 1,4
38 C(K) = .0D0
C
C END CORRECTIONS
NTI = 2*NTI
CRIT = 2.*CRIT
DEN = .0D0
IF (DABS(E2).LT.DE) DEN = DE
CALL SS ( EMIN, DEN, S1, DELTA, EC, BET)
IF (DABS(E2).LT.DE) S12 = S1(2)
CALL TN(1.D0, TRY, EMIN-EF, T)
C(1) = -FMUL*T(1)*S12
DEN = .0D0
IF (DABS(E3).LT.DE) DEN = DE
CALL SS ( EMAX, DEN, S1, DELTA, EC, BET)
IF (DABS(E3).LT.DE) S12 = S1(2)
CALL TN( 1.D0, TRY, EMAX-EF, T)
C(1) = C(1)+FMUL*T(1)*S12
IF (DABS(E2).GT.DE.AND.DABS(E3).GT.DE) GOTO 45
C
C END CORRECTIONS IF E2 OR E3 = 0
CALL TN( 1.D0, TRY, -EF, T)
C(2) = C(2)-PI*FMUL*T(2)/(3.*RP3*B2)
C(1) = C(1)+FMUL*T(1)*PI*(2.*B4/3.D0+.5D0*(RP3+1./B2))/(RP3*EC)
TEMP = -2.D0*FMUL*PI*B4*DE/(3.D0*RP3*EC*DABS(EC))
C(1) = C(1)+TEMP*T(1)
CALL TN (1.D0, TRY, DSIGN(DE,EC)-EF, T)
C(1) = C(1)+4.D0*TEMP*T(1)
CALL TN( 1.D0, TRY, 2.*DSIGN(DE,EC)-EF, T)
C(1) = C(1)+TEMP*T(1)
45 NTI = NTI/2
CRIT = CRIT/2.D0
NL = NT+N-NF
IF (NL.LT.3) GOTO 41
IF (NL.GT.NE) NL = NE
DO 40 L = 1,NL
NTL = NF-N+L
IF (NTL.GT.NT) GOTO 41
TMM = TM(2,NTL)
IF(TRY.NE.0) GOTO 39
IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(2,NTL-NT+5)
39 C(4) = C(4)+S(L,3)*TMM
DO 40 K = 1,2
TMM = TM(K,NTL)
IF (TRY.NE.0) GOTO 40
IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(K,NTL-NT+5)
40 C(K) = C(K)+S(L,K)*TMM
41 DO 42 K = 1,2
42 C(3) = C(3)+C(K)
C(4) = FNE*C(4)
NN=N/NPR
IF (ISIG.GE.1.AND.NN*NPR.EQ.N) print 9800, EF, ( C(K), K = 1,4)
9800 FORMAT (F12.8,4D11.3)
IF (NN*NPR.EQ.N) NAM((LGR-1)*(NF/NPR)+NN) = C(M)
IF (LGR.NE.1.OR.N.NE.1) GO TO 48
CMIN = C(M)
CMAX = C(M)
48 CMIN = DMIN1(CMIN, C(M))
CMAX = DMAX1(CMAX, C(M))
49 CONTINUE
50 EF = EF+DE
C DO 52 K=1,LGR
C NBOUND(K,1) = CMIN
C52 NBOUND(K,2) = CMAX
LGR=LGR+1
IF (ISIG.GE.0) GO TO 22
LGR=LGR-1
CALL GRAF (EFMIN,EFMAX,NAM,NBOUND,NF/NPR,LGR ,FUN)
print 9900, CMIN, CMAX
9900 FORMAT (' CMIN = ', E12.4,' CMAX = ', E12.4)
C
IF (ISIG.EQ.2) GO TO 1
IF (ISIG.EQ.-1) STOP
C
C ERROR MESSAGES
90 print 1100
1100 FORMAT (' NE TOO LARGE')
GOTO 1
93 print 1200
1200 FORMAT(' NT TOO LARGE')
GOTO 22
95 print 1300
1300 FORMAT (' NE TOO SMALL')
GOTO 1
END
C
SUBROUTINE TN(B, TRY, EMU, T )
C
IMPLICIT REAL *8 (A-H,O-Z)
C SUBROUTINE TO CALCULATE THE INTEGRAL OF THE FERMI FUNCTION
C OF E = B*Z**2, T(1), AND ITS DERIVATIVE WITH RESPECT TO
C FERMI LEVEL, T(2), TRY IS THE TEMPERATURE (IN ENERGY UNITS)
C AND EMU = E - MU,WHERE MU IS THE FERMI LEVEL,
DIMENSION T(2)
COMMON/B/ CRIT, DE ,NT
C
PI = 3.1415927D0
IF (DABS(EMU).LT..1*DE) EMU = .0D0
T(1) = .0D0
T(2) = .0D0
N = NT
IF (TRY.EQ.0.) GO TO 40
BETA = 1.D0/TRY
BEMU = BETA*EMU
IF (BEMU.GT.CRIT) GO TO 70
IF (BEMU.LT.-CRIT) GO TO 80
Z0 = 0.D0
ARF = ((-CRIT)/BETA-EMU)/B
C IF (ARF.GT.0D0) Z0 = DSQRT(ARF)
ZL = DSQRT((((CRIT+5.D0)/BETA) -EMU)/B)
C
DZ = (ZL-Z0) /(2*N )
OSC = DZ/3.D0
Z = Z0
N2P1 = 2*N + 1
DO 50 L = 1, N2P1
OSC = -OSC
EX = .0
SIMP = DZ + OSC
IF (L.EQ.1.OR.L.EQ.2*N +1) SIMP = DZ/3.D0
BZ2 = B*Z*Z
ARG = BETA * ( EMU + BZ2 )
IF (DABS(ARG).GT.4.*CRIT) GO TO 50
EX = DEXP(ARG)
D = 1.D0 + EX
DFDM = BETA*EX/D**2
T(1) = T(1) + 2.D0*SIMP*BZ2*DFDM
T(2) = T(2) + SIMP*DFDM
50 Z = Z + DZ
IF (EX.EQ.0D0) RETURN
TEMP = .5/(EX*BETA*B*(Z-DZ))
DO 60 L = 1, 2
T(L) = T(L) + TEMP
60 TEMP = BETA*TEMP
RETURN
C
40 IF (EMU.GE.0.D0) RETURN
T(1) = DSQRT(-EMU/B)
T(2) = -.5D0*T(1)/EMU
RETURN
C
70 IF (BEMU.GT.50) RETURN
EX = DEXP(-BEMU)
TEMP = .5D0*EX*DSQRT(PI/(BETA*B))
DO 75 L = 1, 2
T(L) = TEMP
75 TEMP = TEMP*BETA
RETURN
C
C
80 T(1) =DSQRT(-EMU/B)
T(2) = -.5D0*T(1)/EMU
TEMP = (.5D0*T(2)/EMU)*((PI/BETA)**2)/6.D0
T(1) = T(1) + TEMP
T(2) = T(2) + 1.5*TEMP/EMU
RETURN
END
C
C
SUBROUTINE TONE(X,T)
IMPLICIT REAL *8 (A-H,O-Z)
DIMENSION T(2,4)
X2 = X*X
R = DSQRT(X)
T(1,1) = X2*R*(X2/9.D0-.2D0 )/3.D0
T(2,1) = X*R*(X2/6.D0-X2/7.D0+1.D0/9.D0-1.D0/6.D0)
T(1,2) = X2*R*(-X2/9.D0+X/7.D0+.4D0 )
T(2,2) = X*R*(3.D0*X2/7.D0-2.D0*X/5.D0-2.D0/3.D0-(X+1.)*(X-2.)/2.)
T(1,3) = X*R*(X*X2/9.D0-2.D0*X2/7.D0-X/5.D0+2.D0/3.D0)
T(2,3) = -X*R*(3.D0*X2/7.D0-4.D0*X/5.D0-1.D0/3.D0)
T(1,4) = -X2*R*(X2/9.D0-3.D0*X/7.D0+2.D0/5.D0)/3.D0
T(2,4) = X*R*(3.D0*X2/7.D0-6.D0*X/5.D0+2.D0/3.D0)/3.D0
RETURN
END
SUBROUTINE GRAF (EFMIN,EFMAX,NAM,NBOUND,NF,LGR,FUN)
LOGICAL*1 FIELD,EQSMBL,SMBLS(5)
LOGICAL LEQ,LALLX
INTEGER GR,DIGY,GRINTX,FSTNX
REAL *8 LBX,STX,LBY,SCY,EFMIN,EFMAX
REAL NAM(310),NBOUND(10,2),FUN(NF,LGR)
DATA FIELD/1H /,EQSMBL/1H*/,SMBLS/1H*,1H+,1H=,1H^,1H-/
LALLX=.FALSE.
LEQ=.FALSE.
DO 22 L=1,LGR
DO 22 N=1,NF
22 FUN(N,L)=NAM((L-1)*NF+N)
LBX=EFMIN
STX=(EFMAX-EFMIN)/(NF-1)
LBY=NBOUND(1,1)
SCY=NBOUND(1,2)
CALL MAZIL2(FIELD,EQSMBL,LEQ,LGR,NF,FUN,SMBLS,LBX,STX,
* LBY,SCY,3,3,10,1,LALLX,0)
RETURN
END